!#######################################################
 module iniconds
! Setup of initial conditions for different fields
!#######################################################

!*******************
  contains
!*******************

!========================================
 subroutine set_bc_flow
! Setup velocity boundary conditions
!========================================
   use share; use domain; use flow
   implicit none
   integer :: i
   real (kind=8) :: a
   character (len=100) :: readval
!=========================================

!-----------------
! Top
!-----------------

!  free-slip everywhere
   if(bcT_flow.eq.0) then
     bcT=0

!  constant velocity everywehere
   elseif(bcT_flow.eq.1.and.uTconst.eq.1) then
     bcT=1
     call readinput("top_u",readval)
     read(readval,*) a
     uT=a

!  variable velocity and/or free-slip 
   elseif((bcT_flow.eq.1.and.uTconst.eq.0).or.(bcT_flow.eq.2)) then
     call readinput("top_u_model",readval)
     read(readval,*) i

     if(i.eq.1) then
       call kinematic_slab

     elseif(i.eq.2) then
       call free_slab

     else
       print *, "Top boundary condition model", i, "does not exist. Sorry."
       stop
     endif

   endif  

!-----------------
! Bottom
!-----------------

!  free-slip everywhere
   if(bcB_flow.eq.0) then
     bcB=0

!  constant velocity everywehere
   elseif(bcB_flow.eq.1.and.uBconst.eq.1) then
     bcB=1
     call readinput("bottom_u",readval)
     read(readval,*) a
     uB=a

!  variable velocity and/or free-slip 
   elseif((bcB_flow.eq.1.and.uBconst.eq.0).or.(bcB_flow.eq.2)) then
     print *, "No model for variable bottom boundary velocity has been programmed yet"
     print *, "Cheers"
     stop

   endif  

!-----------------
! Left
!-----------------

!  free-slip everywhere
   if(bcL_flow.eq.0) then
     bcL=0

!  constant velocity everywehere
   elseif(bcL_flow.eq.1.and.vLconst.eq.1) then
     bcL=1
     call readinput("left_v",readval)
     read(readval,*) a
     vL=a

!  variable velocity and/or free-slip 
   elseif((bcL_flow.eq.1.and.vLconst.eq.0).or.(bcL_flow.eq.2)) then
     print *, "No model for variable left boundary velocity has been programmed yet"
     print *, "Cheers"
     stop

   endif  

!-----------------
! Right
!-----------------

!  free-slip everywhere
   if(bcR_flow.eq.0) then
     bcR=0

!  constant velocity everywehere
   elseif(bcR_flow.eq.1.and.vRconst.eq.1) then
     bcR=1
     call readinput("right_v",readval)
     read(readval,*) a
     vR=a

!  variable velocity and/or free-slip 
   elseif((bcR_flow.eq.1.and.vRconst.eq.0).or.(bcR_flow.eq.2)) then
     print *, "No model for variable right boundary velocity has been programmed yet"
     print *, "Cheers"
     stop

   endif  

 end subroutine set_bc_flow


!=========================================================================
 subroutine kinematic_slab
! Setup variable velocity boundary conditions for the top boundary
! NB: the boundary condition array bcT must be set for every i in [1,nx]
!     either to 0 (free-slip) or to 1 (prescribed velocity).
!     At the gridpoints for which bcT(i)=1, the velocity vector must uT(i)
!     must be also set
!=========================================================================
   use share; use domain; use flow
   implicit none
   integer :: i
   real (kind=8) :: u_slab, x_trench
   character (len=100) :: readval
!=========================================
  
   call readinput("u_slab",readval)
   read(readval,*) u_slab

   call readinput("x_trench",readval)
   read(readval,*) x_trench

   do i=1,nx
   
     if(xgrid(i).le.x_trench) then
! Prescribed velocity
       bcT(i)=1
       uT(i)=u_slab

     else
! Prescribed velocity
       bcT(i)=1
       uT(i)=0d0

     endif

   end do

 end subroutine kinematic_slab


!=========================================================================
 subroutine free_slab
! Setup variable velocity boundary conditions for the top boundary
! NB: the boundary condition array bcT must be set for every i in [1,nx]
!     either to 0 (free-slip) or to 1 (prescribed velocity).
!     At the gridpoints for which bcT(i)=1, the velocity vector must uT(i)
!     must be also set
!=========================================================================
   use share; use domain; use flow
   implicit none
   integer :: i
   real (kind=8) :: x_trench
   character (len=100) :: readval
!=========================================

   call readinput("x_trench",readval)
   read(readval,*) x_trench

   do i=1,nx
   
     if(xgrid(i).le.x_trench) then
! Free-slip
       bcT(i)=0

     else
! Prescribed velocity
       bcT(i)=1
       uT(i)=0d0

     endif

   end do
   
 end subroutine free_slab


!===================================
 subroutine initial_temp
! Initial temperature distribution
!===================================
   use share; use domain 
   use markeraux; use flow
   implicit none
   integer :: i,j,wn
   real (kind=8) :: x,z,r,Tini,Tpot,Tadi,top_tbl,bot_tbl,delta,ph,rnd,Tsolidus
   real (kind=8) :: T1,z1,T2,z2,T_bot_tbl,T_top_tbl,Tcool,arg,age_plate,age_over,x_trench
   real (kind=8), dimension(0:3) :: csol
   character (len=100) :: readval
!========================================================================

   call readinput("init_temp_dist",readval)
   read(readval,*) itd   !  initial temperature distribution 

   call readinput("init_temp_pert",readval)
   read(readval,*) delta   !  initial temperature perturbation for sin wave

   call readinput("init_wave_number",readval)
   read(readval,*) wn     !  wave number of the initial sin perturbation

   call readinput("init_wave_phase",readval)
   read(readval,*) ph     !  phase of the initial sin perturbation

   call readinput("temp_rnd_pert",readval)
   read(readval,*) rnd     ! add a small random perturbation 

!----------------------------
! Temperature distributions
!----------------------------

! Linear profile 
   if(itd.eq.1) then

     do i=0,nx+1
       x=xgrid(i)
       do j=0,nz+1
          z=zgrid(j)

           Tnew(i,j)=(1d0-z)
!           Tnew(i,j)=(1d0-z)+delta*dcos(pi/dimx*x)*dsin(pi/dimz*z)

        end do
     end do

! Constant temperature + TBLs 
   elseif(itd.eq.2) then

     call readinput("init_temp",readval)
     read(readval,*) Tini   !  initial internal temperature

     call readinput("top_tbl",readval)
     read(readval,*) top_tbl   ! thickness of the initial top thermal boundary layer

     call readinput("bot_tbl",readval)
     read(readval,*) bot_tbl   ! thickness of the initial bottom thermal boundary layer

     Tnew=0d0
     do j=0,nz+1
       z=zgrid(j)

       if((dimz-z).lt.top_tbl) then
         Tnew(:,j)=Tini*dsin(0.5*pi*(dimz-z)/top_tbl)
       elseif(z.lt.bot_tbl) then
         Tnew(:,j)=1d0-(1d0-Tini)*dsin(0.5*pi*z/bot_tbl)
       else
         Tnew(:,j)=Tini
       endif

     end do

! Solidus + TBLs 
   elseif(itd.eq.3) then

     call readinput("top_tbl",readval)
     read(readval,*) top_tbl   ! thickness of the initial top thermal boundary layer

     call readinput("bot_tbl",readval)
     read(readval,*) bot_tbl   ! thickness of the initial bottom thermal boundary layer

! Non-dimensional coefficients of the solidus
     call readinput("csolidus_0",readval)
     read(readval,*) csol(0)   
     call readinput("csolidus_1",readval)
     read(readval,*) csol(1)   
     call readinput("csolidus_2",readval)
     read(readval,*) csol(2)   
     call readinput("csolidus_3",readval)
     read(readval,*) csol(3)   

     Tnew=0d0
     do j=0,nz+1
       z=zgrid(j)
       
       Tsolidus = 0d0
       do i=0,3
         Tsolidus = Tsolidus + csol(i)*(dimz-z)**dble(i) ! Solidus temperature
       end do

       if((dimz-z).lt.top_tbl) then
         Tnew(:,j)=Tsolidus*dsin(0.5*pi*(dimz-z)/top_tbl)
       elseif(z.lt.bot_tbl) then
         Tnew(:,j)=1d0-(1d0-Tsolidus)*dsin(0.5*pi*z/top_tbl)
       else
         Tnew(:,j)=Tsolidus
       endif

     end do

! Adiabatic profile with a given potential temperature + TBLs 
   elseif(itd.eq.4) then

     call readinput("pot_temp",readval)
     read(readval,*) Tpot   !  initial internal temperature

     call readinput("top_tbl",readval)
     read(readval,*) top_tbl   ! thickness of the initial top thermal boundary layer

     call readinput("bot_tbl",readval)
     read(readval,*) bot_tbl   ! thickness of the initial bottom thermal boundary layer

     Tnew=0d0
     do j=0,nz+1
       z=zgrid(j)

       Tadi=Tsurf*(dexp(Di*(dimz-z))-1d0) + Tpot*(dexp(Di*(dimz-z)))
       
       if((dimz-z).lt.top_tbl) then
         Tnew(:,j)=Tadi*dsin(0.5*pi*(dimz-z)/top_tbl)
       elseif(z.lt.bot_tbl) then
         Tnew(:,j)=1d0-(1d0-Tadi)*dsin(0.5*pi*z/bot_tbl)
       else
         Tnew(:,j)=Tadi
       endif

     end do


! Adiabatic profile with a given potential temperature + half-space cooling (subduction model)
   elseif(itd.eq.5) then
     call readinput("pot_temp",readval)
     read(readval,*) Tpot   !  initial internal temperature
     call readinput("age_over",readval)
     read(readval,*) age_over   !  non-dimensional age of the overriding plate
     call readinput("x_trench",readval)
     read(readval,*) x_trench

     Tnew=0d0
     do j=0,nz+1
      z=zgrid(j)
      Tadi=Tpot*(dexp(Di*(dimz-z)))+Tsurf*(dexp(Di*(dimz-z))-1d0)

      do i=0,nx+1
       x=xgrid(i)
       if(x<x_trench) then
        age_plate=x/x_trench*age_over
       else
        age_plate=age_over
       endif
       arg=(dimz-z)/dsqrt(4*age_plate) ! =(dimz-z)/dsqrt(4*kappa_plate*age_plate), kappa_plate=kappa_0=1d0
       !Tcool=Tsurf + (Tpot-Tsurf)*derf(arg) ! temperature of asthenosphere?
       Tcool=Tpot*derf(arg) ! temperature of asthenosphere?
 
       Tnew(i,j)=Tadi-(Tpot-Tcool)
      end do
     end do

! Temperature prescribed at the base and at the top of the upper and lower TBLs, respectively.
! T increases linearly across both TBLs and across the mantle.
   elseif(itd.eq.6) then

     call readinput("top_tbl",readval)
     read(readval,*) top_tbl   ! thickness of the upper TBL
     call readinput("top_tbl_temp",readval)
     read(readval,*) T_top_tbl   !  temperature at the base of the upper TBL

     call readinput("bot_tbl",readval)
     read(readval,*) bot_tbl   ! thickness of the initial bottom thermal boundary layer
     call readinput("bot_tbl_temp",readval)
     read(readval,*) T_bot_tbl   !  temperature at the top of the lower TBL

     T1 = T_bot_tbl
     z1 = bot_tbl
     T2 = T_top_tbl
     z2 = dimz - top_tbl
     Tnew=0d0
     do j=0,nz+1
       z=zgrid(j)

       if(z.le.z1) then
         Tnew(:,j)= (T1-1d0)/z1*z + 1d0
       elseif(z.ge.z2) then
         Tnew(:,j) = -T2/(1d0-z2)*z + T2/(1d0-z2) 
       else
         Tnew(:,j)= (T1-T2)/(z1-z2)*z + T1 - (T1-T2)/(z1-z2)*z1
       endif

     end do

   endif

!-----------------------
! Perturbations
!-----------------------

! Add a horizontal perturbation
   if(delta.gt.0d0) then
     do i=0,nx+1
       x=xgrid(i)
       do j=0,nz+1
         z=zgrid(j)
         Tnew(i,j)=Tnew(i,j)+delta*dcos(dble(wn)*pi/dimx*x+ph)
       end do
     end do

   endif

! Add a random perturbation
   if(rnd.gt.0d0) then
     r=1d0
     do i=1,nx
       do j=1,nz
         if(mod(j,2).eq.0) then
           call random_number(r)
           Tnew(i,j)=Tnew(i,j)+r*rnd
         else
           call random_number(r)
           Tnew(i,j)=Tnew(i,j)-r*rnd
         endif
       end do
     end do

   endif

! Ensure that BC are satisfied
   Tnew(nx+1,:)=Tnew(nx,:)
   Tnew(0,:)=Tnew(1,:)
   if(bcB_T.eq.0) Tnew(:,0)=1d0
   Tnew(:,nz+1)=0d0

 end subroutine initial_temp

!==================================================
 subroutine initial_melt
! Initialize solidus, liquidus on the grid and
! and depletion both on markes and grid
!==================================================
   use share; use domain 
   use markeraux; use flow
   implicit none
   integer :: k,j
   real (kind=8) :: z
   real (kind=8), dimension(0:3) :: csol,cliq
   character (len=100) :: readval
!========================================================================

!-----------------------------------------------------------------
! Set the initial depletion to 1 both on markers and on the grid
!-----------------------------------------------------------------
   call readinput("depl_mk",readval)
   read(readval,*) depl_mk   ! depletion on markers

   if(depl_mk.eq.1) then 
     deplmk=0d0
     depl=0d0
     depl_old=0d0
   endif

!-----------------------------------
! Set initial solidus and liquidus 
!-----------------------------------
   call readinput("csolidus_0",readval)
   read(readval,*) csol(0)   
   call readinput("csolidus_1",readval)
   read(readval,*) csol(1)   
   call readinput("csolidus_2",readval)
   read(readval,*) csol(2)   
   call readinput("csolidus_3",readval)
   read(readval,*) csol(3)   

   call readinput("cliquidus_0",readval)
   read(readval,*) cliq(0)   
   call readinput("cliquidus_1",readval)
   read(readval,*) cliq(1)   
   call readinput("cliquidus_2",readval)
   read(readval,*) cliq(2)   
   call readinput("cliquidus_3",readval)
   read(readval,*) cliq(3)   

   Tsol=0d0
   Tliq=0d0
   do j=0,nz+1
     z=zgrid(j)
     
     do k=0,3
       Tsol(:,j) = Tsol(:,j) + csol(k)*(dimz-z)**dble(k) ! Solidus
       Tliq(j) = Tliq(j) + cliq(k)*(dimz-z)**dble(k) ! Liquidus
     end do

   end do
   Tsol0(:)=Tsol(1,:)


 end subroutine initial_melt


!==================================================
 subroutine initial_markers
! Initial distribution of markers properties
!==================================================
   use share; use domain 
   use markeraux; use flow
   implicit none
   integer :: i,j,l,nl,k
   integer, dimension(4) :: ai,aj
   real (kind=8) :: xm,zm,dm,r,rnd
   real (kind=8), dimension(4) :: w
   real (kind=8), dimension(0:nx+1,0:nz+1) :: pom
   real (kind=8), dimension(:), allocatable :: zlayer,clayer,hslayer,visclayer
   character (len=1)  :: ch
   character (len=100)  :: readval
   real (kind=8) :: mg_number, delta_rho_gnt_maj, delta_rho_pyx_maj, delta_rho_ol_gamma_ol
   real (kind=8), parameter :: olivine_vol1=0.6d0, pyx_vol1=0.4d0, olivine_vol2=0.56d0, pyx_vol2=0.44d0, garnet_vol=1d0,  &
                               gamma_olivine_vol=0.5d0, majorite_vol=0.5d0, rho_olivine_mg=3227d0, rho_olivine_fe=4402d0, &
                               rho_pyx_mg=3240d0, rho_pyx_fe=3900d0, rho_garnet_mg=3565d0, rho_garnet_fe=4312d0,          &
                               rho_gamma_olivine_mg=3549d0, rho_gamma_olivine_fe=4850d0, rho_majorite_mg=3523d0,          & 
                               rho_majorite_fe=4630d0
!==============================================================================================================================

   call readinput("comp_mk",readval)
   read(readval,*) comp_mk   ! composition on markers

   call readinput("heats_mk",readval)
   read(readval,*) heats_mk   ! heat sources on markers

   call readinput("visc_mk",readval)
   read(readval,*) visc_mk   ! viscosity on markers

   call readinput("markers_interp",readval)
   read(readval,*) mkinterp  ! markers-to-grid weighting scheme

!-----------------------------------------
! Set composition on markers
!-----------------------------------------
   if(comp_mk.eq.1) then 

     if(debug.eq.1) write(*,*) "Set composition on markers"

     call readinput("init_comp_dist",readval)
     read(readval,*) icd   !  initial composition distribution 

     call readinput("comp_rnd_pert",readval)
     read(readval,*) rnd   !  add a small random perturbation 

     Gmk=0d0
!----------------------------------------------------------------------------     
! Linear profile (heavy material on top with G=1, light at the bottom G=0)
!----------------------------------------------------------------------------     
     if(icd.eq.1) then

       do l=1,Nmk
         zm=zmk(l)
         Gmk(l)=zm
       end do

!----------------------------------------------------------------------------     
! Linear profile (heavy material at the bottom with G=1, light at the top G=0)
!----------------------------------------------------------------------------     
      elseif(icd.eq.2) then

       do l=1,Nmk
         zm=dimz-zmk(l)
         Gmk(l)=zm
       end do

!-----------------------------------
! Different compositional layers
!-----------------------------------
     elseif(icd.eq.3) then

       call readinput("n_comp_layers",readval)
       read(readval,*) nl !  number of different compositional layers
       allocate(zlayer(nl+1),clayer(nl))
       zlayer(nl+1)=dimz

       do k=1,nl
         write (ch, '(I1)') k
         call readinput("clayer_z_"//ch,readval)
         read(readval,*) zlayer(k)
         call readinput("clayer_C_"//ch,readval)
         read(readval,*) clayer(k)
       end do

       do l=1,Nmk
         zm=zmk(l)

         do k=1,nl
           if(zm.ge.zlayer(k).and.zm.le.zlayer(k+1)) then
             Gmk(l)=clayer(k)
             mktype(l)=k
           endif
         end do

       end do
       deallocate(zlayer,clayer)

!---------------------------------------------------------------------
! Magma ocean overturn profile (without phase transition at 14 GPa)
! see Elkins-Tanton et al., Meteoritics and planetary science (2003)
!---------------------------------------------------------------------
     elseif(icd.eq.4) then
       do l=1,Nmk
         dm=dimz-zmk(l) !depth

         if(dm.le.1d0.and.dm.gt.0.578d0) then
           Gmk(l) = -0.27891d0*dm + 0.3371d0
         endif

         if(dm.le.0.578d0.and.dm.gt.0.528d0) then
           Gmk(l) = 0.3824d0
         endif

         if(dm.le.0.528d0.and.dm.gt.0.133d0) then
           Gmk(l) =  0.377829d0 - 1.19154*dm + 0.901422*dm*dm
         endif

         if(dm.le.0.133d0.and.dm.gt.0d0) then 
           Gmk(l) =  1d0 - 10.3443d0*dm + 34.5463*dm*dm
         endif

       end do

!----------------------------------------------------------------------
! Magma ocean overturn profile (with phase transition at 14 GPa)
! see Elkins-Tanton et al., Meteoritics and planetary science (2003)
!----------------------------------------------------------------------
     elseif(icd.eq.5) then

       do i=1,Nmk
         dm=dimz-zmk(i) !depth

         if((dm.le.1d0).and.(dm.gt.0.578)) then 
               mg_number = (76.89099 + 7.109*dm)/100.0
               delta_rho_gnt_maj = (mg_number*rho_majorite_mg + (1.0-mg_number)*rho_majorite_fe) - &
                                   (mg_number*rho_garnet_mg + (1.0-mg_number)*rho_garnet_fe)
               delta_rho_pyx_maj = (mg_number*rho_majorite_mg + (1.0-mg_number)*rho_majorite_fe) - &
                                   (mg_number*rho_pyx_mg + (1.0-mg_number)*rho_pyx_fe) 
               delta_rho_ol_gamma_ol = (mg_number*rho_gamma_olivine_mg + (1.0-mg_number)*rho_gamma_olivine_fe) - &
                                       (mg_number*rho_olivine_mg + (1.0-mg_number)*rho_olivine_fe)

               Gmk(i) =  gamma_olivine_vol*mg_number*rho_gamma_olivine_mg +                                       &
                         gamma_olivine_vol*(1.0-mg_number)*rho_gamma_olivine_fe +                                 &
                         majorite_vol*mg_number*rho_majorite_mg + majorite_vol*(1.0-mg_number)*rho_majorite_fe
               Gmk(i) = (Gmk(i) - 3330.0)/620.0

               dGmk(i) = (0.2*delta_rho_gnt_maj + 0.3*delta_rho_pyx_maj + 0.5*delta_rho_ol_gamma_ol)
               dGmk(i) = dGmk(i)/620.0

               mktype(i) = 2
         endif

         if ((dm.le.0.578).and.(dm.gt.0.528)) then 
               mg_number = 0.85
               delta_rho_gnt_maj = (mg_number*rho_majorite_mg + (1.0-mg_number)*rho_majorite_fe) - &
                                   (mg_number*rho_garnet_mg + (1.0-mg_number)*rho_garnet_fe)
               delta_rho_pyx_maj = (mg_number*rho_majorite_mg + (1.0-mg_number)*rho_majorite_fe) - &
                                   (mg_number*rho_pyx_mg + (1.0-mg_number)*rho_pyx_fe) 
               delta_rho_ol_gamma_ol = (mg_number*rho_gamma_olivine_mg + (1.0-mg_number)*rho_gamma_olivine_fe) - &
                                       (mg_number*rho_olivine_mg + (1.0-mg_number)*rho_olivine_fe)

               Gmk(i) = garnet_vol*mg_number*rho_garnet_mg + garnet_vol*(1.0-mg_number)*rho_garnet_fe 
               Gmk(i) = (Gmk(i) - 3330.0)/620.0

               dGmk(i) = delta_rho_gnt_maj
               dGmk(i) = dGmk(i)/620.0

               mktype(i) = 1
         endif

         if ((dm.le.0.528).and.(dm.gt.0.133)) then 
               mg_number = (83.65316 + 10.12658*dm)/100.0
               delta_rho_gnt_maj = (mg_number*rho_majorite_mg + (1.0-mg_number)*rho_majorite_fe) - &
                                   (mg_number*rho_garnet_mg + (1.0-mg_number)*rho_garnet_fe)
               delta_rho_pyx_maj = (mg_number*rho_majorite_mg + (1.0-mg_number)*rho_majorite_fe) - &
                                   (mg_number*rho_pyx_mg + (1.0-mg_number)*rho_pyx_fe) 
               delta_rho_ol_gamma_ol = (mg_number*rho_gamma_olivine_mg + (1.0-mg_number)*rho_gamma_olivine_fe) - &
                                       (mg_number*rho_olivine_mg + (1.0-mg_number)*rho_olivine_fe)

               Gmk(i) = olivine_vol2*mg_number*rho_olivine_mg + olivine_vol2*(1.0-mg_number)*rho_olivine_fe + &
                        pyx_vol2*mg_number*rho_pyx_mg + pyx_vol2*(1.0-mg_number)*rho_pyx_fe
               Gmk(i) = (Gmk(i) - 3330.0)/620.0

               dGmk(i) = 0.44*delta_rho_pyx_maj + 0.56*delta_rho_ol_gamma_ol
               dGmk(i) = dGmk(i)/620.0

               mktype(i) = 1
         endif

         if ((dm.le.0.133).and.(dm.ge.0.0)) then 
               mg_number = (25.0 + 451.1278*dm)/100.0
               delta_rho_gnt_maj = (mg_number*rho_majorite_mg + (1.0-mg_number)*rho_majorite_fe) - &
                                   (mg_number*rho_garnet_mg + (1.0-mg_number)*rho_garnet_fe)
               delta_rho_pyx_maj = (mg_number*rho_majorite_mg + (1.0-mg_number)*rho_majorite_fe) - &
                                   (mg_number*rho_pyx_mg + (1.0-mg_number)*rho_pyx_fe) 
               delta_rho_ol_gamma_ol = (mg_number*rho_gamma_olivine_mg + (1.0-mg_number)*rho_gamma_olivine_fe) - &
                                       (mg_number*rho_olivine_mg + (1.0-mg_number)*rho_olivine_fe)

               Gmk(i) = olivine_vol1*mg_number*rho_olivine_mg + olivine_vol1*(1.0-mg_number)*rho_olivine_fe + &
                        pyx_vol1*mg_number*rho_pyx_mg + pyx_vol1*(1.0-mg_number)*rho_pyx_fe
               Gmk(i) = (Gmk(i) - 3330.0)/620.0
               
               dGmk(i) = 0.4*delta_rho_pyx_maj + 0.6*delta_rho_ol_gamma_ol
               dGmk(i) = dGmk(i)/620.0

               mktype(i) = 1
         endif

       end do

     endif  

   endif

!-----------------------------------------
! Set RaH on markers
! NB: It overwrites RaH read from yacc.inp
!-----------------------------------------
   if(heats_mk.eq.1) then 
     RaHmk=0d0

     if(debug.eq.1) write(*,*) "Set RaH on markers"

     call readinput("init_heat_dist",readval)
     read(readval,*) ihd   !  initial distribution of heat sources

! Linear distribution growing from the bottom with height (RaH=0 at the bottom,
! RaH=RaH at the top)
     if(ihd.eq.1) then

       do l=1,Nmk
         zm=zmk(l)
         RaHmk(l)=RaH(1,1)*zm
       end do

! Distribution in different layers
     elseif(ihd.eq.2) then
       call readinput("n_heats_layers",readval)
       read(readval,*) nl !  number of different layers

       allocate(zlayer(nl+1),hslayer(nl))
       zlayer(nl+1)=dimz

       do k=1,nl
         write (ch, '(I1)') k
         call readinput("heatslayer_z_"//ch,readval)
         read(readval,*) zlayer(k)
         call readinput("heatslayer_RaH_"//ch,readval)
         read(readval,*) hslayer(k)
       end do

       do l=1,Nmk
         zm=zmk(l)
  
         do k=1,nl
           if(zm.ge.zlayer(k).and.zm.le.zlayer(k+1)) then
             RaHmk(l)=hslayer(k)
             mktype(l)=k
           endif
         end do

       end do

       deallocate(zlayer,hslayer)
     endif

   endif

!-----------------------------------------------------------
! Set viscosity eta on markers
! NB: It overwrites eta calculated from subroutine varparm
!-----------------------------------------------------------
   if(visc_mk.eq.1) then 
     etamk=0d0

     if(debug.eq.1) write(*,*) "Set viscosity on markers"

     call readinput("init_visc_dist",readval)
     read(readval,*) ivd   !  initial distribution of viscosity

! Distribution in different layers
     if(ivd.eq.1) then
       call readinput("n_visc_layers",readval)
       read(readval,*) nl !  number of different layers

       allocate(zlayer(nl+1),visclayer(nl))
       zlayer(nl+1)=dimz

       do k=1,nl
         write (ch, '(I1)') k
         call readinput("visclayer_z_"//ch,readval)
         read(readval,*) zlayer(k)
         call readinput("visclayer_eta_"//ch,readval)
         read(readval,*) visclayer(k)
       end do

       do l=1,Nmk
         zm=zmk(l)
  
         do k=1,nl
           if(zm.ge.zlayer(k).and.zm.le.zlayer(k+1)) then
             etamk(l)=visclayer(k)
           endif
         end do

       end do

       deallocate(zlayer,visclayer)
     endif

   endif

!---------------------------------------------------------------------
! Calculate initial concentration 'Gam' by interpolating the marker 
! concentration onto the grid
!---------------------------------------------------------------------
   if(comp_mk.eq.1) Gam=0d0
   if(heats_mk.eq.1) RaH=0d0
   if(visc_mk.eq.1) eta=0d0

   pom=0d0
   do l=1,Nmk
     xm=xmk(l)
     zm=zmk(l)

!  Compute the weights of all markers found in the 4 cells
!  that surround each grid point 
    call Wmarker4(xm,zm,ai,aj,w)
    do i=1,4 
      pom(ai(i),aj(i))=pom(ai(i),aj(i))+w(i)
      if(comp_mk.eq.1) Gam(ai(i),aj(i))=Gam(ai(i),aj(i))+Gmk(l)*w(i)
      if(heats_mk.eq.1) RaH(ai(i),aj(i))=RaH(ai(i),aj(i))+RaHmk(l)*w(i)
      if(visc_mk.eq.1) eta(ai(i),aj(i))=eta(ai(i),aj(i))+etamk(l)*w(i)
    end do  

   end do
   if(comp_mk.eq.1) Gam=Gam/pom
   if(heats_mk.eq.1) RaH=RaH/pom
   if(visc_mk.eq.1) eta=eta/pom
    
   if(rnd.gt.0d0) then
! Add a small random perturbation to the composition field
     r=1d0
     do i=1,nx
       do j=1,nz
         if(mod(j,2).eq.0) then
           call random_number(r)
           Gam(i,j)=Gam(i,j)+r*rnd
         else
           call random_number(r)
           Gam(i,j)=Gam(i,j)-r*rnd
         endif
       end do
     end do
   endif

   if(comp_mk.eq.1)  Gam_old=Gam
   if(heats_mk.eq.1) RaH_old=RaH
   if(visc_mk.eq.1)  eta_old=eta

 end subroutine initial_markers


!#######################################
 end module iniconds
!####################################### 

